library(Zelig)
library("ZeligChoice")
library(foreign)
library(xtable)

#Cleaning up data 

im<-read.dta("data_final.dta",convert.dates=T)
im$imop<-as.factor(im$imop)
im$ctry<-as.factor(im$ctry)
im$year<-as.factor(im$year)
im.simp<-im[which(im$hiallsample==1 &  im$native==1 & im$retired==0),]

imclean<-subset(im.simp,select=c("imop",  "imchg","empchg", "isced3",
                     "isced4","isced5","essround",
                     "age","male","union","ess2Xproximity2", 
                     "ess3Xproximity2", "ess4Xproximity2",
                     "ctry","year","proximity2","pcforborn","urate","benefits",
                     "gdppercapitaconstant2000us","imoe"))
imclean<-imclean[which(imclean$year!=2010),]
imclean<-imclean[which(imclean$year!=2011),]
imclean<-imclean[which(imclean$year!=2012),]
imclean<-imclean[which(imclean$year!=0),]
imclean<-na.omit(imclean)
imclean$year<-factor(imclean$year)
imclean$ctry<-factor(imclean$ctry)
imclean$year<-factor(imclean$year)
imclean$essround<-factor(imclean$essround)
table(imclean$year)


#For all countries

z.out<- zelig(imop ~ ess4Xproximity2+empchg+ isced3 +isced4 + isced5+
                age+male+union+essround+ess2Xproximity2+ 
                ess3Xproximity2+
                ctry+proximity2+urate+gdppercapitaconstant2000us+
                pcforborn+benefits,
              model ="oprobit", data = imclean)
summary(z.out)
z.out$vcov

###########################################################################
#Plotting
##########################################################################
library(gplots)

#Figure 1

x.low<-setx(z.out,ess4Xproximity2=0)

preds<-0
my.preds<-c()
my.sds<-c()
exp<-c()
sds<-c()

for(i in 1:60){
  preds<-preds+.1
  x.high<-setx(z.out,ess4Xproximity2=preds)
  s.out <- sim(z.out,x=x.low,x1=x.high)
  exp[i]<-s.out$stats$`Expected Values: P(Y=j|X1)`[1,1] 
  sds[i]<-s.out$stats$`Expected Values: P(Y=j|X1)`[1,2] 
}

a<-seq(from=.1,to=6,by=.1)
a<-exp(a)

par(mfrow=c(1,2))
plotCI(x=exp[1:40], y = a[1:40], uiw=sds[1:40],
       col="darkblue",barcol="blue",pch=16,gap=0,sfrac=0,
       main="Elections and 'Allow None'",
       xlab="Expected Value of 'allow none'",
       ylab="Days to Election",bty="n",err="x")

plot(exp[1:35],a[1:35],pch=16,col="darkblue",main="30 Days to Election",
     xlab="Expected Value of 'allow none'",
     ylab="days to election",bty="n")

#Figure 2

x.low<-setx(z.out,ess4Xproximity2=0,essround=4)
x.high<-setx(z.out,ess4Xproximity2=6,essround=4)
s.out <- sim(z.out,x=x.low,x1=x.high)
par(mfrow=c(1,1))
plot(s.out)
s.out$stat

#Figure 3

x.low<-setx(z.out,essround=2,ess4Xproximity2=1)
x.high<-setx(z.out,essround=4,ess4Xproximity2=1)
s.out <- sim(z.out,x=x.low,x1=x.high)
par(mfrow=c(1,1))
plot(s.out)
s.out$stat

#Figure 4

z.out<- zelig(imoe ~ isced3 +isced4 + isced5+
                age+male+union+essround+ess2Xproximity2+ 
                ess3Xproximity2+ ess4Xproximity2+proximity2+
                ctry,
              model ="oprobit", data = imclean)
summary(z.out)


x.low<-setx(z.out,ess4Xproximity2=0)

preds<-0
my.preds<-c()
my.sds<-c()
exp<-c()
sds<-c()

for(i in 1:60){
  preds<-preds+.1
  x.high<-setx(z.out,ess4Xproximity2=preds)
  s.out <- sim(z.out,x=x.low,x1=x.high)
  exp[i]<-s.out$stats$`Expected Values: P(Y=j|X1)`[1,1] 
  sds[i]<-s.out$stats$`Expected Values: P(Y=j|X1)`[1,2] 
}

a<-seq(from=.1,to=6,by=.1)
a<-exp(a)

par(mfrow=c(1,2))
plotCI(x=exp, y = a, uiw=sds,
       col="darkblue",barcol="blue",pch=16,gap=0,sfrac=0,
       main="Elections and 'Allow None' of Same Ethnic Group",
       xlab="Expected Value of 'allow none'",
       ylab="days to election",bty="n",err="x")

plot(exp[1:35],a[1:35],pch=16,col="darkblue",main="30 Days to Election",
     xlab="Expected Value of 'allow none'",
     ylab="days to election",bty="n")

#figure 5

imclean8<-imclean[which(imclean$ctry=="IE"),]
imclean8$essround<-factor(imclean8$essround)

z8<- zelig(imoe ~ proximity2,
           model ="oprobit", data = imclean8[which(imclean8$essround==3),])
summary(z8)
x.low<-setx(z8,proximity2=0)
x.high<-setx(z8,proximity2=3)
s.out <- sim(z8,x=x.low,x1=x.high)
summary(s.out)
plot(s.out)


preds<-0
my.preds<-c()
my.sds<-c()
exp<-c()
sds<-c()

for(i in 1:60){
  preds<-preds+.1
  x.high<-setx(z8,proximity2=preds)
  s.out <- sim(z8,x=x.low,x1=x.high)
  exp[i]<-s.out$stats$`Expected Values: P(Y=j|X1)`[1,1] 
  sds[i]<-s.out$stats$`Expected Values: P(Y=j|X1)`[1,2] 
}
a<-seq(from=.1,to=6,by=.1)
a<-exp(a)
par(mfrow=c(1,1))
plotCI(x=exp, y = a, uiw=sds,
       col="darkgreen",barcol="green",pch=16,gap=0,sfrac=0,
       main="Elections and 'Allow None' in Ireland",
       xlab="Expected Values of 'allow none'",
       ylab="Days to Election",bty="n",err="x")
plotCI(x=exp[1:50], y = a[1:50], uiw=sds,
       col="darkblue",barcol="blue",pch=16,gap=0,sfrac=0,
       main="Elections and 'Allow None'",
       xlab="Expected Value of 'allow none'",
       ylab="days to election",bty="n",err="x")
